## Replication of the Democratic governed state analysis from:
##
## Christopher Adolph, Kenya Amano, Bree Bang-Jensen, 
## Nancy Fullman, Beatrice Magistro, Grace Reinke, Rachel Castellano,## Megan Erickson, and John Wilkerson, "The Pandemic Policy U-Turn:## The role of partisanship, public health, and race in decisions
## to ease COVID-19 social distancing policies in the U.S.", 
## forthcoming in Perspectives on Politics.
##
## This file:
## make Figure 4 and computes the contents of Table B2
##
## Estimate pooled stratified Cox proportional hazards of 
## the first indoor easing of state-level social distancing measures, 
## 16 April to 6 July 2020; and produce plots of hazard rates 
## and expected acceleration FOR DEMOCRATIC GOVERNED STATES ONLY
##
## Christopher Adolph
## 18 June 2021     faculty.washington.edu/cadolph

## Clear memory
rm(list=ls())

## Set random number seed for consistency in replication
set.seed(123456)

## Use fancy fonts?
fancy <- FALSE

## Load libraries
library(survival)      # For Cox ph and simulation
library(coxed)         # for AME predicted survival times
## note that a helper file is loaded to fix a typo in
## coxed2; see below
library(MASS)          # for mvrnorm()
library(formula.tools) # for lhs(), rhs()
library(tile)          # graphics; available from
## http://faculty.washington.edu/cadolph/software
library(RColorBrewer)  # Color selection
library(simcf)         # for extractdata; available from
## http://faculty.washington.edu/cadolph/software

## Load helper functions, including a fix for an apparent bug in coxed 0.3.3
##  (see helperAnalysis.R for more details)
source("helperAnalysis.R")

## Set up colors
col <- brewer.pal(7, "Set1")
nicered <- col[1]
niceblue <- col[2]
nicegreen <- col[3]
nicepurple <- col[4]
niceorange <- col[5]
nicebrown <- col[7]

## Load data
load("AnalysisData.Rdata", verbose=TRUE)

## Save new copy of data for additional processing for the baseline model
allsurvData <- mData

## Keep only states with Democratic governors
## This is the key difference from the baseline analysis
allsurvData <- allsurvData[allsurvData$DemGovernor==1,]

## Set earliest date = April 16 >=50 (check)
startDate <- 50
allsurvData <- allsurvData[(allsurvData$DaysSinceFCT>=startDate), ]

## Set end date = July 6 <=131 (check)
endDate <- 131
allsurvData <- allsurvData[(allsurvData$DaysSinceFCT<=endDate), ]

## Set source of COVID case, death, and test data
allsurvData <- allsurvData[allsurvData$COVIDSource=="JHU",]  # was JHU

## Select policy types
allsurvData <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                           (allsurvData$PolicyType=="AnyBusinessClose")|
                           (allsurvData$PolicyType=="RestaurantRestrict")|
                           (allsurvData$PolicyType=="BarRestrict")|
                           (allsurvData$PolicyType=="StayAtHome")
                          ,]

## Set up model specification, with clustering and stratification
model1a <- sdEasing ~
    TrumpVote2016 +
    PctBlack +
    TestPositivityRateLagSlopeLM +
    logNewDeathsPerM7DaysMAVG +
    NoDeaths7DaysMAVG +
    NewConfirmedCasesPerMLagSlopeLM +
    log(PopDensity) +
    strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) +
    cluster(as.factor(StateName))

## Create a dataset with only at risk cases
obsSurvData <- allsurvData[!is.na(allsurvData$PolicySubNonReligIndoorEaseDV),]

## Create survival object
sdEasing <-  Surv(time=obsSurvData$DaysSinceFCT-1,
                  time2=obsSurvData$DaysSinceFCT,
                  event=obsSurvData$PolicySubNonReligIndoorEaseDV
)

## Renumber rows (required by coxed below)
rownames(obsSurvData) <- 1:nrow(obsSurvData)

## Estimate CPH model: this is the baseline model
res1a <- coxph(model1a, data=obsSurvData, x=TRUE, y=TRUE)

## This summary is used to fill in the goodness of fit for Table B2.
## However, the counterfactuals for covariates come from the
## simulations below for Figure 4, as they involve continuous
## covariates, for which the one-unit-change hazard ratios are unhelpful
summary(res1a)
AIC(res1a)

###################################################################################
## Set up counterfactuals for Figure 4

## Coefficient numbers
xNum <- c(1, 2, 3, 4, 6, 7)

## Scenario names
scenname <- c("More Trump Voters",
              "Higher Black Population",
              "Downward Trend Positivity",
              "Lower Rate of New Deaths",
              "Downward Trend New Cases",
              "Lower Population Density"
              )

## Note: the next two commands must deal with an issue in computing quantiles 
## (in particular, the 25th percentile) created by replacing 0s with 1s 
## in the unlogged death rate data:  We use the uncorrected origNewDeathsPerM7DaysMAVG
## to compute the correct quantiles, while preserve the replaced values in 
## logNewDeathsPerM7DaysMAVG for in-sample simulation via coxed
            
## "Before" scenario (theoretically discourages easing)
## Unlogged versions of the below are also used 
## to fill in the "pre" values in Table B2
xpre1 <- with(obsSurvData,
               c(quantile(TrumpVote2016, probs=0.25),
                 quantile(PctBlack, probs=0.25),
                 quantile(TestPositivityRateLagSlopeLM, probs=0.75),
                 log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.75)),
                 quantile(NewConfirmedCasesPerMLagSlopeLM, probs=0.75),
                 log(quantile(PopDensity, probs=0.75))
                 ))                

## "After" scenario (theoretically encourages easing)
## Unlogged versions of the below are also used 
## to fill in the "post" values in Table B2
xpost1 <- with(obsSurvData,
              c(quantile(TrumpVote2016, probs=0.75),
                quantile(PctBlack, probs=0.75),
                quantile(TestPositivityRateLagSlopeLM, probs=0.25),
                log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.25)),
                quantile(NewConfirmedCasesPerMLagSlopeLM, probs=0.25),
                log(quantile(PopDensity, probs=0.25))
                ))

## Make Figure 4; which must then be polished in Adobe Illustrator.
## The object hr is also used to fill in Table B2 hazard rates 
## and confidence intervals, and to label points in the polished plots
hr <- plotCovidCPH("Figure4",  res1=res1a,
                   scenname=scenname, xNum=xNum, xpre1=xpre1, xpost1=xpost1, 
                   limits=c(0.1, 3.1),
                   at=c(0.5, 1, 1.5, 2, 2.5,3.0),
                   labsX=c("0.50", "1.00","1.50", "2.00", "2.50", "3.00"),
                   labsT=list(c("half", "as", "likely"),
                              c("just", "as", "likely"),
                              c("50%", "more", "likely"),
                              c("100%", "more", "likely"),
                              c("150%", "more", "likely"),
                              c("200%", "more", "likely")),
                   sortBy=1, rev=TRUE, logaxis=FALSE,
                   titleT="The first move to ease social distancing is...",
                   titleX="Hazard Ratio",
                   pch=c(16, 16, 15, 16, 17, 16),
                   col=c(niceblue, "black", nicered, niceblue, nicepurple, niceblue),
                   fancy=fancy
                   )
